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PREFACE 


EUs  Memorandum  documents^ a  numerical  procedure  uhich  vas 

* 

demised  to  solve  the  problem  of  hypervelocity  impact  yid  Vhieh  was  # 

* 

subsequently  revised  for  application  to  the  problem  of  wrateriftg  and 

Ground  shock  from  a  nuclear  surface  burst. 

The  RAND  symposium  on  Hich-Speed  Impact,  held  in  1355#  focused 

attention  on  a  forthcoming  Air  Force  need  for  infonnation  relative 

to  hypervelocity  impact.  In  fact,  this  symposium  va s  the  first  of 

a  series  of  six  symposia  sponsored  Jointly  by  the  Air  Force,  Army,  . 

•  • 

and  Navy  on  the  sane  topic  since  that  date. 

At  the  first  symposium,  technical  data  vere  presented  vhieh 
strongly  suggested  that  the  hypervelocity-impaot  process  vas  hydro- 
dynamic  fa  nature,  involving  the  substantial  compression  of  even  the 
strongest  materials,  leading  to  shooks  and  severe  fluid  distortion 
in  the  resulting  flow.  Hie  problem  vas  hopelessly  complicated  from 
•  the  analytical  point  of  view,  and  at  that  time,  numerical  techniques 
did  not  exist  which  were  adequate  to  provide  the  desired  solutions* 
The  numerical  procedure  discussed  in  this  Iteaoranttom  vas  devised 
specifically  to  attack  this  problem. 

On§e  developed,  the  method  proved  to  have  application  to  an 
area  of  Alar  Force  interest  not  eor  template!  in  the  original  research, 
namely  era  taring  and  ground  dhoel;  induced  by  a  nuclear  surface  hurst* 
Ohe  procedure  has  been  applied  to  this  prcWen,  and  the  results  care 
documented  in  RM-fcobO,  Crater  inn  From  ^teteten  Surface  &tcst* 


o  © 
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3UMMARY 


This  Memorandum  discusses  In  detail  a  numerical  method  for  sclv- 

•  • 

ing  the  compressible,  hydrodynamic  equations  under  the* limitations 

•  • 

of  (1)  two  space  dimension,  (2)  the  inviscid  approximation,  and  (3) 

the  adiabatic  approximation.  The  method  allows  for  the  occurrence 

» 

of  shocks,  contact  discontinuities,  and  interfaces.  Under  a  proper 

•  • 

prescription  of  initial  and  boundary  conditions,  the  method  generates 

solutions  including  the  above  physical  phenomena. 

♦ 

The  basis  of  the  method  is  the  extension  to  two  space  dimensions 

* 

of  the  particle- in- cell  (PIC)  concept  first  proposed  by  Harlow  for  a  ♦ 
one -dimensional  computational  scheme.  The  PIC  concept  involves  mass 
points  moving  in  a  Lagrangian  sense  through  an  Euler i an  space  grid, 
Besides  mass,  the  points  carry  with  them  the  proper  amount  of  momentum, 
kinetic  energy,  and  internal  energy. 

Ifce  computational  method  approximates  a  set  of  partial  differen¬ 
tial  equations  containing  terms  in  addition  to  those  of  the  compressible, 
hydrodynamic  equations  under  the  approximations  cited.  The  terms  are 
qualitatively  suggestive  of  thermal  conductivity  and  viscosity  but 
are  l*!>t  exactly  analogous  to  them.  Hie  terras  are  responsible  for  the 
Computational  stability  of  the  numerical  method^  since  they  smear 
ghock  fronts  over  a  few  grid  spaces.  * 

the  evaluation  of  errors  due  to  these  terms  ftnd  tljuse  of  higher  • 

Order  has  not  a£  yet  yielded  to  analysis.  The  errors  have  been 

♦ 

* 

eYaltteoed  by  comparing  numerical  solution#  with  analytical  solutions 

Of  test  problems  and  with  numerical  solutions  of  one -dimensional 

•  •  •  . 

problems .  •  • 
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The  computational  scheme  contains  a  feature  known  as  "grid- 
changing,”  which  permits  optimum  resolution  of  all  hhases  of  the  prob¬ 
lem  using  the  limited  memory  capacity  of  present-day  electronic  com¬ 
puters  * 

•  •  • 
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e  ■  specific  internal  energy 

g  ■  gravitational  acceleration  0 
P  -  pressure  # 
t  -  time 

u  *  radial  component  of  particle  velocity 

u  *  particle  vector  velocity  •  • 

e 

V  -  volume  • 

v  m  axi al#c omponant  of  particle  velocity 

vg  *  velocity  of  sound 
x  -  radial  coordinate 
y  ■  axial  coordinate 
cp  -  stress  power 

4 

$  *  potential  of  the  external  force  field 

1  -* 

t  =  total  specific  energy,  e  +  j  u  *  u 

•  *  i 
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i,  nfTROptygrion 

In  recent  times,  many  important  physical  problems  that  demand 
the  solution  of  the  compressible  hydrodynamic  equations,  using  more 
than  me  space  dimension,  have  arisen.  Two  examples  are  the  ground 
motion  in  the  early  stages  of  a  nuclear  ground  burst  and  the  motion 
Induced  when  a  projectile  strikes  a  target  at  extremely  high  veloci¬ 
ties,  It  might  seem  curious  that  although  both  problems  Involve 
solid  media,  the  hydrodynamic  equations  are  used  to  describe  the 
phenomena  However,  this  is  precisely  the  approximation  found 
necessary  in  the  high-pressure,  high -density  regions  encountered 
In  these  processes. 

The  solutions  of  this  class  of  physical  problems  generally 
feature  the  presence  of  shocks,  which  are  allowed  In  the  framework 
of  the  nonlinear  partial  differential  equations.  Across  the "e 
shocks,  discontinuities  in  the  dependent  variables  occur.  The  com¬ 
plications  which  arise  In  the  analytical  treatment  of  shooks  are 
generally  very  great.  For  this  reason  a  great  deal  of  effeyt  has 
been  devoted  in  the  past  to  the  formulation  of  schemes  that  yield 
numerical  solutions  on  electronic  computers.  Because  of  these  com¬ 
plications,  most  of  the  work  in  the  past  has  been  confined  to  prob¬ 
lems  which  contain  a  single  space  variable.  This  in  turn  Implies 
the  existence  of  a  strong  spatial  symmetry,  usually  plane,  spherical, 
*  or  cylindrical,  • 

However,  in  the  case  of  solids,  the  presence  of  a  free  surface 

often  has  an  extremely  important  Influence  on  physical* behavior. 

Such  a  free  surface  amplifies  the  complexity  of  the  already  difficult 
•  * 


e 
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analytleal  treatment  to  such  a  degree  that  workers  In  the  field 

despair  of  attaining  analytical  solutions  by  the  use  of  presently 

known  techniques,  A  free  surfaoe  also  Increases  the  complexity  of 

the  numerical  solutions  In  that  It  demands  the  use  of  more  than  one 
* 

space  dimension,  and  Its  presence  visually  leads  to  severe  distortion 
If  the  medium, 

«  # 

At  the  time  that  this  work  was  undertaken,  it  was  evident  that 
the  numerical  analogue  of  the  Lag rang i an  formulation  broke  down  when 
severe  distortion  of  the  medium  occurred  and  that  the  analogue  of 
the  Euler lan  formulation  contained  spurious  material-diffusion  terms 
whenever  the  problem  included  free  surfaces  and  interfaces.  The 
technique  described  in  this  Memorandum  circumvents  these  two  diffi¬ 
culties,  Furthermore,  it  has  been  successfully  applied  to  the 

hypervelocity- Impact  problem^  and  to  the  early  motion  of  the  ground 

•  (2) 

•  during  a  nuclear  surface  burst,  ' 

The  example  chosen  far  describing  the  method  is  that  of  the  air 

0 

flow  following  a  nuclear  airburst.  It  should  be  emphasized  that 
this  problem  has  not  yet  been  solved  by  the  method,  and  that  diffi¬ 
culties  that  prevent  its  solution  may  arise.  However,  this  problem 
contains  some  generalizations  not  required  in  the  problems  that  have 
been  solved,  thereby  allowing  a  more  complete  discussion  to  be  given. 
These  are  the  presence  of  an  external  body  force  (gravity)  and  a 
nonhcnogeneous  medium  (the  exponentially  varying  atmosphere). 
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II .  PBOBUM  DESCRIPTIOH 

SYSTEM  OF  EQUATIONS  TO  BE  SOLVED 

®ie  numerical  method  is  designed  to  solve  the  hydrodynamic 
equations  with  the  inviscid,  adiabatic  approximation.  The  represen¬ 
tation  of  these  equations  in  Eulerian  coordinates  takes  the  follow- 
lag  form: 


pjj-  +  grad  P  +  p  grad  4*0  0) 

p  div  ^  ^ 

2®  ?  2e  =,  o  (3) 

Dt  ^5  Dt 

P  »  P(p,  e)  ^ 


vhere 


0 

Dt 


b 

& 


+  u  •  grad 


The  Independent  variables  are  the  time  t,  and  a  set  of  spatial  coor 
dinates  x.  The  dependent  variables  are 


u  *  particle  velocity 
p  «  pressure 
p  «=  density 

e  *  specific  internal  energy 


The  potential  *  of  the  external  force  field  must  be  specified  in  advance. 
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Equation  (l),  Euler's  equation  of  motion,  contains  the  assvap- 
tion  that  the  only  forces  which  accelerate  the  fluid  are  pressure 
forces  and  external  forces  idiich  nay  he  derived  from  a  potential. 

Equation  (2),  the  equation  of  continuity.  Is  a  mathematical 
statement  of  the  fact  that  mass  sust  he  conserved. 

Equation  (3)  Is  the  first  lav  of  themodynaad.es  under  the  adia¬ 
batic  approximation.  It  states  that  the  only  vay  the  Internal  energy 
of  a  fluid  element  may  be  changed  is  through  the  action  of  pressure 
forces  during  expansion  or  congress  ion  of  the  element. 

Equation  (4)  Is  the  equation  of  state  of  the  substance  under 
consideration.  It  establishes  an  equilibrium  relation  (solved 
explicitly  far  P)  between  the  pressure,  density,  and  specific  Inter¬ 
nal  energy  of  a  small  element  of  the  material. 

Although  Eqs.  (l)  through  (4),  together  with  a  potential  func¬ 
tion  #  and  appropriate  boundary  and  Initial  conditions,  completely 
specify  the  motion,  they  are  not  In  the  most  convenient  form  for 
the  particular  m*erical  computations  vhlch  we  have  In  mind.  Ve 
therefore  transform  them  In  the  following  vay. 

Hie  dot  product  of  Eq.  (l)  with  u  is  taken  to  obtain 

pu  •  ^  +  u  •  grad  P  +  pu  •  grad  4=0  (5) 

Substitution  of  Eq.  (2)  into  Eq.  (3)  yields 

-e§S  “  p  div  u  (6) 


We  then  use  the  vector  identity 


0 
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dlv(Pu)  -  P  dlv  u  +  u  •  grad  P 
together  with  Eq.  (6),  to  obtain 

u  *  grad  P  «  div(Pu)  +  pg|  (7) 

Equation  (7)  is  substituted  into  Eq.  (?)•  The  result  is 

pj-g  |^u  •  u  +  e  1  +  u1y(f5)  +  pu  •  grad  *  -  0  (8) 

If  we  integrate  Eq.  (8)  over  a  region  whose  volume  is  V  and  whose 
surface  is  denoted  by  S,  we  have 

^  [p3  •  grad  *  +  u  •  u  +  e^J  dV  •  -J  Pu  •  d3  (9) 

where  Gauss's  theorem  has  been  applied  to  the  right-hand  side. 

Equation  (9)  Is  the  form  which  is  actually  used  in  the  numeri¬ 
cal  computations,  together  with  Eqs.  (l)  and  (U). 

DflTIAL  AND  BOUNDARY  CONDITIOHS 

Ike  problems  which  we  wish  to  solve  are  transient,  nonlinear, 
initial-boundary  value  problems  in  two  space  dimensions.  Their 
nature  is  such  that  if  we  prescribe  initial  conditions  at  time  t 
and.  appropriate  boundary  conditions  for  all  times  T,  where 

t  s  t  1  t  +  At 

we  can  find  solution  fleldr  at  a  subsequent  time  t  +  At.  As  we 
have  stated  earlier,  such  problems  do  not  lend  themselves  to  an 
analytical  treatment,  especially  when  the  solution  fields  contain 
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a  variety  of  discontinuities  such  as  shocks,  slip  lines,  and  free 

•ft 

surfaces,  as  they  do  in  many  problems  of  Interest* 

We  are  therefore  forced  to  solve  difference  analogues  of  the 
differential  equations  to  obtain  solutions  at  successive  points  In 
time* 

In  order  to  solve  the  difference  equations  at  a  succession  of 
times  it  Is  necessary  to  have  at  our  disposal  the  following  kinds 
of  data  t&ich  must  be  used  simultaneously  with  the  difference  equa¬ 
tions* 

Initial  Data 

At  some  time  tQ  values  of  all  of  the  dependent  variables,  (P, 
p,  e,  u)  must  be  specified  at  a  set  of  points  covering  the  region 
of  Interest.  It  is  not  necessary  to  specify  all  these  quantities  at 
the  same  set  of  points*  In  the  numerical  scheme  which  Is  descr  bed 
later,  the  velocity,  density,  and  specific  Internal  energy  are  speci¬ 
fied  at  one  set  of  points,  while  the  pressure  Is  specified  at  an 
entirely  different  set  of  points*  The  reason  for  this  will  became 
man  evident  later* 

Boundary  Data 

Boundary-data  'requirements  assise  a  variety  of  forms,  i&lch 
depend  on  the  kinds  of  boundarlec  that  appear  In  the  region  we  wish 
to  cover  by  our  mserlcal  solution* 


- V - 

At  the  present  time  we  have  not  even  the  assurance  from  mathe¬ 
matics  that  these  problems  are  solvable  and  contain  unique  solu¬ 
tions.^  Sene  assurance  is  obtained  by  cooparlx^  numerical  solutions 
with  physical  experiments,  with  other  types  of  numerical  solutions, 
and  with  the  very  few  analytical  solutions  which  can  be  obtained. 
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First,  there  art  rigid  boundaries.  The  appropriate  constraint 
htrt  It  that  tht  normal  o  mponent  of  the  fluid  velocity  vast  coin¬ 
cide  with  the  normal  component  of  velocity  assigned  to  tut  rigid 
boundary.  If  the  rigid  boundary  la  fixed  in  tlms,  than  the  ocm&l 
component  of  fluid  Telocity  must  vanish,  of  course,  at  this  boundary. 
If  a  point,  axis,  or  plane  of  i/nwtry  appears  in  the  problem,  It 
must  be  treated  as  a  fixed  boundary. 

A  second  kind  of  boundary  often  encountered  is  a  free  surface. 

At  a  free  surface  we  must  assign  zero  pressure  to  the  fluid,  other- 
vise  It  vould  do  vork  on  a  vac  mm  or  vice  versa. 

At  Interfaces  and  contact  discontinuities  the  pressure  sad  mor¬ 
ns!  components  of  fluid  velocity  must  be  continuous,  vhile  density, 
internal  energy,  and.  tangential  components  of  velocity  may  change 
abruptly  across  such  surface j  •  In  the  mmsrlcal  procedure  to  be 
described,  no  special  arrangements  need  to  be  made  in  order  to  pre¬ 
serve  the  requisite  continuities  In  pressure  and  normal  velocities. 
The  principal  reason  far  this  is  that  the  numerical  procedure  intro¬ 
duces  terms  vhlch  have  the  effect  of  smoothing  out  discontinuities 
(k  5) 

In  the  solution  fields.  *  However,  because  we  keep  track  of 
each  particle  In  the  system,  we  can  still  trace  out  a  fairly  sharp 
"Interface"  vhlch  separates  wave  elements  comprised  of  different 
materials. 

Contact  discontinuities  can  be  approximately  determined  by 
noticing  the  distribution  of  mass  points  for  a  single  material  at 
a  given  tima.^ 

Finally,  moving  shock  surfaces  nay  occur  on  vhlch  the  well- 
known  Ranklne-Hugonlot  conditions  must  be  s&tlsflel.  Again,  err  eg  a 


o 
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in  the  numerical  pirooedure  tend  to  replace  a  shock  discontinuity  by 
a  zone  of  finite  width  across  which  the  dependent  variables  vary 
in  a  rapid,  but  relatively  smooth,  fashion.  3%e  Rankine-Hugoniot 
conditions  are  automatically  satisfied  by  the  difference  equations 
governing  the  motion,  so  that  again  no  special  arranges nt  must  be 
made  for  shocks  in  the  numerical  procedure. 
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as  Internal  and  kinetic  energv,  Since  the  particles  move  in  a 
Lagrangian  sense  throughout  an  Eulerian  grid,  the  method  may  he 
said  to  incorporate  some  features  of  both  formulations*  The  depen¬ 
dent  variables,  viz*,  density,  velocity,  pressure,  and  specific 
internal  energy,  are  always  specified  at  the  Eulerian  grid  points. 

In  solving  the  equations,  initial  and  necessary  boundary  con¬ 
ditions  are  set  up.  The  differential  equations  are  then  integrated 
with  respect  to  time .  Thus  the  essential  problem  is,  given  the 
variable  fields  at  time  t,  to  find  the  new  fields  at  time  t  +  fit.  The 
computation  proceeds  in  two  steps.  In  the  first,  the  convective 
terms  (starred  in  Eqs*  (10a)  through  (10c))  are  neglected,  and  the 
truncated  set  of  differential  equations  is  approximated  in  finite 
difference  form.  From  these  it  is  possible  to  obtain  new  tentative 
values  of  velocity  and  specific  internal  energy,  u  and  e*,  respec¬ 
tively* 

In  the  second  step,  the  mass  points  are  moved,  using  the  average 
of  the  old  and  the  tentative  new  velocity.  After  the  mass  movement, 
a  calculation  is  made  to  determine  which  masses  have  changed  cells. 

If  no  masses  enter  a  cell  during  this  time  cycle,  the  tentative 
values  are  accepted  as  the  final  values. 

When  one  or  more  masses  enter  a  cell,  a  process  known  as 
repartitioning  is  carried  out.  Tne  particle  is  first  considered  to 
have  brought  with  it  an  amount  of  momentum,  given  by  the  product 
of  its  mass  and  the  velocity  of  the  cell  which  it  left.  This  momen¬ 
tum  is  added  to  that  of  the  cell  the  mass  entered,  and  a  new  velocity 
for  that  cell  is  determined  by  dividing  the  total  new  momentum  by 
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the  total  new  mass.  The  new  velocity  Is  taken  to  be  the  final 
velocity  for  this  particular  tine  cycle;  and  the  nev  total  mass,  the 
final  mass.  In  addition,  a  mass  entering  a  cell  Is  assxned  to  have 
brought  with  It  an  amount  of  Interna],  energy,  given,  similarly,  by 
the  product  of  Its  mass  and  the  specific  Internal  energy  of  the  cell 
vhich  It  left.  This  Is  added  to  the  Internal  energy  of  the  cell 
which  the  mass  entered. 

It  Is  readily  shovn  that  the  repartition! ng  process  does  not 
conserve  kinetic  energy  If  the  velocity  of  the  cell  which  the  particle 
entered  ic  different  from  that  which  It  left.  When  the  moment  ia  and 
interned  energy  are  conserved  In  this  way,  same  kinetic  energy  Is 
always  lost.  That  Is,  the  estimated  kinetic  energy  assigned  to  the 
cell  after  mass  movement  but  before  repartition! ng  Is  always  greater 
than  the  kinetic  energy  of  the  cell  after  repartitioning .  In  order 
to  conserve  total  energy,  the  loss  of  kinetic  energy  Is  arbitrarily 
added  to  the  Internal  energy  of  the  cell  In  which  this  kinetic- energy 
defect  occurred.  The  mass  movement  and  subsequent  repartltlonlng 
process  simulates  the  convective  terns  nagleeted  In  the  first  step. 

DISCU5SI01!  OF  ERRORS 

A  thorough  discussion  of  FTC  errors  would  include  the  determi¬ 
nation  of  a  single  set  of  difference  and  differential  equations  which 
correspond  to  the  FTC  process,  the  establishment  of  the  convergence 
of  PIC  solutions  to  solutions  of  the  fluid-dynamic  equations  and 
their  Initial  and  boundary  data,  and  the  establishment  of  a  stabil¬ 
ity  criterion. 


o 
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At  the  present  time,  there  appears  to  be  no  satisfactory  analy¬ 
tical  treatment  of  the  PIC  convergence  and  stability  characteristics. 
The  greatest  effort  has  been  placed  on  the  determination  of  a  set 
of  difference  equations  which  form  an  analogue  to  the  PIC  numerical 
scheme  and  on  the  construction  of  differential  equations  that  are 
modeled  by  this  process  to  first  order  in  the  time  increment  and 
to  second  order  in  the  cell  widths. 

The  derivation  of  the  differential  equations  which  are  approxi¬ 
mately  represented  by  PIC  solutions  lies  outside  of  the  scope  of 
this  Memorandum.  Those  who  are  interested  in  this  development  may 
refer  to  Refs.  4  and  8. 

For  the  present  numerical  scheme,  these  differential  equations 
have  the  follcving  form  in  cylindrical  coordinates  with  the  angular 
dependence  omitted: 


H  *  i  ^  ^  ' 


(continuity ) 


+  3x  +  ^  =  x  ^  (xXu  $  +  ay  (Xv  $  (momentum,  x  direction) 

+  iy  v  3  *4  (xX«  +  4  (Xv  §>  (momentun,>  y  Erection) 

Do  fl  dxu  &vl  1  d  /  dex  b,  dex 

pDt  p  |_x  Sx  iyj  ”  x  &x  *^u  dx'  dy  '*v  dy' 

-  \  [<S>s  *  <I>S]*  v  [<!>' 


+  l(l^)  +  (^)  (energy) 


vhert. 
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"  P|ul  * 

\  -  pm5! 

In  general;  boundary. cell  analysis  does  not  yield  analogues  to 
these  differential  equations*  Ve  see  that  the  error  terms,  that  1b 
the  terms  which  do  not  coincide  with  our  original  equations  of 
motion;  formally  appear  In  the  above  set  of  equations  In  the  same 
manner  as  vould  the  effects  of  certain  physical  dissipative  mecha¬ 
nisms.  This  has  led  to  the  common  usage  of  the  terminology  "effec¬ 
tive  viscosity"  and  "effective  heat  conduction"  to  describe  these 
error* ,  Moreover,  the  effects  of  these  errors  on  the  computer  solu¬ 
tions  have  been  Interpreted  as  those  which  would.  have  been  caused 
by  an  effective  viscosity  axd  heat  conduction  Introduced  Into  the 
difference  analogues  of  the  original  differential  equations. Ii 
the  right-hand  side  of  the  analogue  moment  us  equation  is  to  be  the 
divergence  of  a  collection  of  elements,  t1^  ,  and  the  terms  on  the 
right-hand  side  of  the  analogue  energy  equation  not  containing  e  are 

- V - 

In  this  section,  liberal  use  of  the  tensor  notation  viil  be 
11  p 

made.  Thus  t  represents  a  collection  of  n  elements  which  trans¬ 
form  like  the  component s  of  a  contravarlant  tensor  of  rank  two;  u1 
and  u^  represent  the  contra  variant  and  covariant  elements,  respec¬ 
tively,  of  the  velocity  vector.  A  coma  Is  used  to  denote  covariant 
differentiation,  and  sun&tlon  Is  assused  (unless  explicitly  stated 

otherwise)  whenever  an  index  Is  repeated,  for  Instance,  9  *  u.  . 

n  .  .  du ,  n  "I  * " 

is  a  shorthand,  far  9  »  X  \  -E  r?.u  where  is  the 

i,J«l  ol  iJ  aJ 

appropriate  Christoff  el  symbol. 


to  be  the  stress  power,  <?,  associated  vlth  t 


U 


*  ■  'lN.j 

where  the  for*  of  neither  of  these  relations  Is  to  explicitly  depend 
on  the  details  of  the  Initial  and  boundary  data  of  the  flow  problem, 
then  these  i»rtlcular  error  te «s  In  the  analogue  equations  wy  be 
associated  viota  a  particular  dynamic- stress  matrix  t(1J)  glTen  by 


<w>  -  3D  *•*  '  -  2  “  *  3)* 

T(3k)  -  t(k3)  •  0  it  -  1,  8,  3 


vhere 


and 


p  I  I  doc^ 

Xj  -  e— 5 - 

{u1,  MZ,  U3}  -  {u,  V,  0^ 
{x\  X2,  X3}  -  (x,  y,  e} 


- V  the  transformation  properties  of  a  set  of  elements  are  not 

known,  th^do  not  satisfy  a  tensor  transformation  law,  we  use 

the  notation  A(lj)  rather  than  AlJ,  A^,  etc.  If  a  repeated  Index 

1.  not  «med  we  put  parenthesis  around  It  to  emphasize  this  fact. 

l^ot  sunned;  however,  A<1J)  ^  represents 

a  sun. 
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Trom  this  foam  far  the  elenente  of  the  eetrlx  y(1J),  end  from  their 
Interpretation  as  « lease ts  of  a  tyaealo- stress  Matrix,  the  follow- 
lag  facts  My  be  1 —ill Italy  derived: 

The  element#  of  t(1J)  do  not  transform  Ilka  the  elements  of  a 
oontravariant  tensor  of  rank  two.  Thus,  the  stress  power  t(lj)  u^j 
is  not  invariant  under  admissible  coordinate  transformations  in 
Thli  leads  to  the  result  that  suoh  physical  properties  as  entropy, 
speclfio  internal  energy,  and  density  at  a  point  will  vary  with  a 
chaise  of  spatial  coordinates  describing  that  points  Moreover,  the 
stress  vector,  F*(p),  at  a  given  point  (p)  with  respect  to  &  surface 
through  the  point  described  by  its  unit  normal  vector  Vj(p),  pre¬ 
sumably  la  related  to  r(ij)  by  F*(p)  -  t(ij)  Vj(p).  Again,  the  lack 
of  tensorlal  property  of  t(ij)  leads  to  F**s  which  depend  not  only 
on  the  surface  v1(p)  but  on  the  particular  coordinate  system  used 
to  describe  the  position  of  points  in  the  portion  of  under  con¬ 
sideration.  (This  is  not  surprising,  since  the  entire  erxer  analy¬ 
sis  leading  to  the  analogue  differential  equation  is  valid  for  one 
particular  coordinate  system  and  no  other.) 

The  established  relations  for  t(1,J)  do  not  leave  the  flow 

« 

equations  Galilean  invariant. 

By  inspection,  we  see  that  t(1J)  Is  not  syamatric  with  respect 
to  interchange  of  the  indices  i  and  J.  If  we  perform  the  usual 
decomposites  into  symmetric  (s)  and  antisymmetric  (A)  parts,  we 
obtain  t(1J)  -  t3(1J)  ♦  ta(1J).  We  know  that  a  constitutive  relation 

■■  M 

This  particular  property  of  the  analogue  equations  has  already 

(5) 

been  recognised  by  Harlow.  ' 
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c  or  responding  to  the  viscosity  hypothesis  vill  only  be  related  to 
tg(ij).  The  antisymmetric  part  T^(i j)  must  be  accounted  far  in 
some  other  way. 

The  viscosity  hypothesis  is  tq(1J)  ■  C(ljkl)  «.  where 
T  1  3  ^ 

€k/  *5  |_\l  +  Ujekl  azki  tiie  C'8  ixr*  symmetrical  vith  respect  to 

interchange  of  the  first  two  and  last  two  indices.  Because  Tg(iJ) 

i  Ikl 

is  not  a  tensor,  ve  avoid  the  usual  representation,  C  ,  which 
implies  tensorlal  properties  that  our  coefficients  do  not  have..  It 
is  easily  shown  that  for  the  given  form  of  Tq(ij)  the  viscosity 
coefficient  C(l,  2,  1,  2)  (and  all  of  those  derivable  from  it  by 
symmetry  arguments)  must  satisfy  a  set  of  equations  which  are  incon¬ 
sistent  except  on  a  line  given  by  J  j  Ax'*'  -  |  u2  |  Ax2 . 

Hence,  the  origin  of  Ts(ij)  cannot  be  said  to  be  viscous  in 
nature.  If  we  introduce  the  most  general  linear  constitutive  rela¬ 
tion  connecting  the  3tress  elements  T ( ij )  vith  the  deformation  ten¬ 
sor  «^,  T(ij)  =  B(ijfci)  we  still  obtain  incanpatiole  relations 

for  some  of  the  coefficients  of  this  linear  relation.  Thus  the 
antisymmetric  part  of  t(1J)  (as  well  as  the  off-diagonal  terms  of 
the  syncetrlc  part)  cannot  be  obtained  by  even  the  most  general 
linear  constitutive  relation  connecting  the  stress  elements  with 
the  deformation  tensor  c 

The  presence  of  the  antisymmetric  elements  T^(ij)  adds  addi¬ 
tional  confusion  to  the  attempted  representation  of  errors  as  physi¬ 
cal  effects,  for  it  can  be  shown  that  a  rather  arbitrary  partition¬ 
ing  among  a  number  of  physical  effects  (couple- stresses,  couple 
density  fields,  a*jd  possibly  more  general  constitutive  relations 


than  we  have  considered  here)  may  all  lead  to  the  same  TA(tj)  that 
is  implied  by  our  analogue  differential  equations.  Thus,  among  other 
difficulties,  there  seems  to  he  no  unique  way  of  assigning  physical 
effects  to  the  error  terms  considered, 

A  similar  analysis  can  he  made  which  start*  with  the  assumption 
that  the  remaining  error  term  in  the  energy  equation  containing  the 
specific  internal  energy,  e,  may  he  thought  of  as  the  divergence  of 


a  heat- flux  vector 


This  analysis  also  exhibits  the  many  weak¬ 


nesses  of  such  an  assumption.  We  first  observe  that  the  elements  of 

the  ’’pseudo-heat- flux  vector”  Q(i),  which  are  given  by  i )}  = 

f\  — —  o,  ol,  do  not  transform  like  the  components  of  a  contra- 

\  1  fcc  d  cbc  J 

variant  vector,  with  the  result  that  the  heat  flux,  as  well  as  the 
physical  variables,  entropy,  specific  internal  energy,  density,  etc., 
exhibit  the  highly  undesirable  property  of  changing  at  a  point  under 
admissible  coordinate  transformations.  Another  drawback  of  this  assump¬ 
tion  is  that  it  implies  a  he  at -conduct  ion  law  where  heat  flow  is  pro¬ 
portional  to  the  gradient  of  specific  internal  energy  rather  than  of 
temperature,  a  condition  which  would  make  physical  senue  only  when  e 
was  Just  a  function  of  temperature  (which  in  general  is  not  true)* 

In  summary,  we  may  say  that  the  assumption  that  the  enors  may 
be  thought  of  as  related  to  distinct  physical  effects  leads  to  many 
contradictions  of  the  assumption.  Therefore  there  seems  to  be 
little  merit  in  treating  the  lowest- order  errors  iaiuced  by  a  PIC 
program  as  anything  other  than  Just  plain  errors. 

Although  the  terms  which  we  have  been  discussing  represent 
errors,  they  are  largely  responsible  for  the  applicability  of  PIC 
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as  a  practical  numerical  scheme.  These  elements  provide  the 
required  smoothing  that  ve  need  if  we  do  not  wish  to  keep  track  of 
shocks  which  would  be  expected  to  appear  in  the  actual  solution  of 
the  inviscid  equations  for  adiabatic  compressive  flows /9>10) 

As  a  result  of  these  terms,  our  solutions  exhibit  relatively 
smooth  variations  in  shocked  regions,  rather  than  the  violent 
oscillations  which  would  characterize  numerical  solutions  in  these 
same  regions  when  no  such  elements  are  present.  These  artificial 
smoothing  errors  are  much  larger  than  real  dissipative  terms  which 
may  occur  in  our  physical  system.  (They  must  be  in  order  to  provide 
the  necessary  stabilizing  influence.)  For  this  reason,  certain  pre¬ 
cautions  must  be  observed  when  using  this  sort  of  a  program.  Shocks 
are  smeared  out  over  approximately  three  grid  increments  in  the 
numerical  solutions,  whereas  the  true  shock  thickness  may  occur 
ever  a  small  fraction  of  a  grid  increment.  Hence,  shock  thicknesses 
derived  from  this  sort  of  program  should  not  be  quoted  as  physi¬ 
cally  significant  results.  Also,  real  dissipative  phenomena  should 
not  be  incorporated  into  the  original  equations  where  the  corre spool¬ 
ing  artificial  smoothing  effects  are  many  times  larger  than  the  real 
ones. 

In  many  instances,  it  is  possible  to  replace  these  artificial 
smoothing  terms  by  more  desirable  forms.  Such  replacements  are 
treated  in  Ref.  5* 

It  is  easy  to  show  that  whenever  mass  points  leave  a  cell  but 


none  enter  it,  the  difference  equations  for  that  cell  contain  neither 
our  previously  defined  errors  nor  any  convective  terms. ^ 
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Thls  particular  situation  seldom  causes  any  poxbles  In  Interior 
cells,  but  It  often  arises  and  persists  for  many  time  cycles  In 
cells  adjoining  a  boundary,  particularly  a  free  surface  in  a  region 
of  essentially  one -dimensional  flow*  The  effects  which  this  situa¬ 
tion  has  on  the  PIC  solution  of  problems  are  little  understood.  It 
suffices  to  say  here  that  such  cells  may  contribute  a  major  soiree 
of  errors  in  PIC  solutions  and  that  these  errors  may  propagate  to 
the  interior  cells* 

Because  of  the  lack  of  a  good  error  analysis  for  this  kind  of 
numerical  scheme,  it  must  be  emphasized  that  the  principal  evaluation 
of  the  results  must  be  based  on  a  comparison  of  PIC  solutions,  and 
mathematical  solutions  whenever  they  are  available  (which  is  seldom). 
Other  tests,  such  as  checking  the  sphericity  of  results  (where  spher¬ 
icity  is  expected)  and  varying  time  and  space  increments  in  the  PIC 
program  to  see  what  happens,  are  sometimes  useful. 

An  empirical  investigation  of  this  kind  was  used  to  ascertain 
the  nature  and  extent  of  errors  incurred  by  the  PIC  numerical  scheme 
such  as  describe!!  in  this  HemorandUD,  and  the  following  tests  were 
made, 

a.  Numerical  solutions  to  one-dimensional,  uniform- shock-wave 
problems  were  solve.’  fear  a  variety  of  shock  strengths  and  for  equa¬ 
tions  of  state  of  alminum,  iron,  tuff,  and  polytropic  gases  with 
gaama  equal  to  1*4  and  2,0*  Numerical  results  were  compared  to  the 
simple  analytical  solutions  which  in  these  instances  were  at  our 
disposal* 
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b.  Numerical  solutions  to  the  problems  of  the  one -dimensional 
rarefaction  of  air  irto  a  vacuum  and  of  compressed  tuff  into  a 
vacuum  were  obtained.  In  each  instance  the  resulting  theoretical 
flows  are  simple  centered  waves,  which  admit  analytical  solutions. 
Again  the  numerical  results  were  compared  to  these  analytical  solu¬ 
tions. 

c.  An  axisymsetric  numerical  solution  to  a  spherical  air-burst 
problem  in  a  uniform  atmosphere  was  computed  in  cylindrical  coordi¬ 
nates  and  compared  to  a  Lagrangian  solution  of  the  same  problem  by 
H.  L.  Brode^10^  in  spherical  coordinates. 

The  results  of  these  tests  led  to  the  following  set  of  general 
qualitative  observations  which  appeared  to  be  consistent  with  our 
numerical  solutions. 

Shock  velocities  and  the  estimated  position  of  shock  fronts 
were  in  good  agreement  with  the  corresponding  theoretical  solutions 
and  with  values  obtained  in  Ref.  10.  Generally,  agreement  within 
±  5  per  cent  was  obtained.  The  estimated  position  of  shock  fronts 
in  the  numerical  solutions  was  taken  to  be  the  point  at  which  the 
pressure  was  equal  to  the  arithmetical  average  of  the  peak  pressure 
behind  the  shock  zone  and  the  initial  pressure.  The  width  of  the 
shock  zones  varied  from  two  to  five  cell  widths  but  generally 
remained  around  three  cell  widths  in  extent. 

The  largest  errors  in  the  solution  fields  generally  occurred 
immediately  behind  a  shock  zone  and  in  regions  of  very  low  density 
where  there  were  very  few  (generally  less  than  four)  masr  points 
per  cell.  Most  often  (though  this  was  not  always  true)  the  numerical 
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solutlon  tended  to  exceed  the  corresponding  theoretical  values  right 
behind  the  shock  and  then  to  execute  a  damped  oscillatory  action 
about  the  theoretical  values.  The  velocity  and  specific  Internal 
energy  generally  exhibited  smaller  errors  than  did  the  density  and 
pressure  in  these  regions.  Relative  errors  seldom  exceeded  20 
per  cent  and  vere  more  often  on  the  order  of  5  per  cent  for 
reasonable  choices  of  the  cell  widths  and  time  increments*  In 
the  continuous  flows  (rarefaction  solutions)  the  numerical  solv 
tions  were  very  good,  often  giving  relative  errors  of  less  than 
1  per  cent  and  exhibiting  stronger  deviations  only  when  regions  of 
very  low  density  were  encountered  in  the  tail  of  the  rarefaction 
for  air.  (The  tuff  rarefaction  exhibited  no  tailing  off  of  the 
density  to  zero.) 

GRID  CHANGES 

Far  such  purely  practical  reasons  as  the  limited  size  of  the 
electronic  computer's  fast  memory  and  the  desire  to  keep  the  comput¬ 
ing  time  within  reasonable  bounds,  it  is  possible  to  use  only  a 
limited  number  of  grid  points.  Consequently,  resolution  in  the  solu¬ 
tion  of  two-dimensional  problems  is  very  troublesome. 

To  utilize  the  resolution  to  the  utmost,  an  artifice  termed  a 
"grid  change”  has  been  developed.  Ve  take  advantage  of  the  fact 
that  the  interesting  effects  expand  in  size  and  that  the  medium 
outside  the  motion  is  undisturbed  until  seme  signal  from  the  expand¬ 
ing  region  reaches  it.  In  the  case  of  hypervelocity  Impact  and 
nuclear  bursts,  the  exterior  fluid  is  undisturbed  until  a  shock 
reaches  it,  and  the  region  of  motion  is  entirely  within  this 


outward-moving  shock.  Thus,  It  Is  only  necessary  that  the  grid 
cover  the  region  of  notion  and  a  reasonable  amount  of  undisturbed 
fluid  outside  of  it.  The  procedure  Is  to  lay  down  such  a  6rld  and 
constantly  test  for  the  slightest  indication  of  notion  at  the  grid's 
p eri^^ry.  The  computation  Is  stopped  when  such  motion  Is  detected, 
and  a  nev,  somewhat  larger  grid  is  laid  dovn  for  the  next  jiiase  of 
the  computation,  -which  encompasses  same  nev  undisturoed  fluid.  In 
this  way,  the  interesting  region  always  occupies  a  large  fraction 
of  the  grid  points  available. 
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IV,  DETAILED  DESQRIPTIOH  OF  THE  METHOD 
PHYSICAL  EXAMPLE 

In  order  to  illustrate  the  numerical  process,  ve  take  as  an 
example  a  surface  explosion  In  air  over  a  fixed  ground  in  the  pre¬ 
sence  of  a  uniform  gravitational  field.  The  geometric  representa¬ 
tion  selected  is  one  of  cylindrical  symmetry  about  the  axis  x  •  0 
with  the  ground  at  y  -  0  as  shown  in  Fig,  1.  Ve  treat  the  air  as 
a  perfect  gas,  so  that  the  equation  of  state  for  the  fluid  Is 

P  (p,  «)  -  (Y  -  1)P«  (11) 

where  y  is  a  constant  usually  equal  to  1.4*  It  will  be  apparent 
that  any  equation  of  state  of  the  form  P  «  f (p,  e)  may  as  readily 
be  used.  The  external  force  field  has  a  potential 

*  ■  gy  (12) 

where  g  is  the  gravitational  acceleration. 

The  initial  conditions  are 

u(x,  o)  «  0  (13) 

p(x,  o)  -  p0  exp(-y/X)  (ik) 

e(x,  o)  -  aeQ  (for  x2  +  y2  *  r2) 
e(x,  o)  -  eQ 


(for  x2  +  y2  >  r2) 


(15) 
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Fig.  1  —  Surface-burst  geometry 


-2?. 


The  particle  Telocity  u  Is  xero  at  t  »  0;  the  density  p  of  the  fluid 
la  an  exponential  function  of  altitude;  the  specific  internal  energy 
e  Is  Cq  In  the  undisturbed  area,  and  a  constant  a  tines  e^  vlthln 
the  disturbed  area.  The  Initial  pressure  is  determined  by  substi¬ 
tution  of  Eqs.  (14)  and  (15)  into  Eq.  (ll),  the  equation  of  state. 

The  boundary  conditions  to  be  satisfied  are  that  there  is  to 
be  no  fluid  motion  across  the  fixed  boundary  y  -  0  and  that  the 
radial  component  of  reloclty  u  vanishes  at  x  ■  0.  These  tvo  con* 
dltlons  may  be  written  as 

y(x,  o,  t)  -  0 
«(o,  y,  t)  -  0 

where  ▼  and  u  are  respectively  the  axial  and  radial  components  of 
velocity.  The  fluid  is  considered  to  extend  indefinitely  with 
increasing  x  and  lncreaslx^  y. 

At  this  point  it  should  be  emphasized  that  the  particular  repre¬ 
sentation  of  the  physical  problem  described  above  Is  not  necessarily 
the  most  realistic  model  for  an  extensive  investigation  of  the  physi¬ 
cal  phenomena;  rather,  it  is  selected  to  illustrate  the  application 
of  the  numerical  procedure.  It  is  clear  that  other  Initial  condi¬ 
tions  vlthln  the  disturbed  region,  other  equations  of  state,  other 
models  of  the  undisturbed  medium,  and  other  heights  of  tarst  may  be 
used  In  place  of  those  assumed  above. 

RgjgjgggmOg  OF  THE  FLUID 

The  volume  Is  considered  to  be  divided  Into  a  finite  number  of 
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volune  elements,  each  of  which  is  fixed  in  position,  shape,  and  size* 
For  our  example  this  region  is  hounded  by  x  *  0  on  the  left,  y  -  0 
on  the  lover  side,  and  extends  indefinitely  in  the  direction  of 
increasing  x  and  y*  Since  ve  wish  to  maintain  a  oertain  degree  of 
resolution,  for  the  time  being,  ve  limit  the  region  to  a  maximum 
on  the  right  and  y^^  above,  as  shown  in  Fig*  1. 

One  of  the  volume  elements  into  which  the  region  is  subdivided 
is  shown  in  Big.  2a.  Element  ij  is  bounded  on  the  left  by  x^,  on 
the  right  by  on  the  lower  side  by  y^,  and  on  the  upper  side 

by  y1J+1*  It  subtends  an  angle  of  one  radian  about  x  *  0*  Ihe 
volume  of  element  ij  is 

vij  “  2  ^xn+i  "  xn^  (yij+i  '  yij) 

A  convenient  notation  is  that  the  center  of  area  of  the  volume  ele¬ 
ment  cross  section  shown  in  Fig*  2b  is  at  x2i,  y^  where 

X21  “  2  (xll  +  Xli+l) 

y2J  -  2  (ylJ  +  ylJ+l) 
and  the  center  of  volume  is  at  x^,  y^  where 


y3J  “  i  (yij  +  yij+i) 


Fig.  2a — Volume  element  Ij  In  cylindrical  symmetry 


Furthermore,  the  subcelle  01 J,  11  J,  21 J,  and  31 J  all  have  equal 
voltaet. 

The  fluid  In  each  volume  element,  or  cell.  Is  represented  by  a 
number  of  discrete  mass  points*  These  mass  points  are  free  to  move 
from  cell  to  cell  during  the  history  of  the  motion.  At  any  instant 
in  time  tQ,  the  total  mass  in  a  given  cell  is  equal  to  the  sum 
of  the  values  of  the  discrete  masses  contained  in  the  cell.  With 
the  total  mass  and  the  volume  of  a  cell  known,  the  density  p,  or  mass 
per  unit  volume,  in  cell  ij  is 


*ij  -  «?Aj 


vhere  is  the  density  at  time  tQ,  aul  is  given  by  Eq.  (16)* 

Velocity,,  density,  and  specific  internal  energy  are  properties 
attached  to  the  cell,  and  they  are  considered  to  be  constant  over 
the  cell  at  a  particular  time.  The  components  of  velocity  at  time 
tQ  are  designated  in  the  x  direction  by  ujj^,  and  in  the  y  direction 
by  v£j;  the  specific  internal  energy  e  for  cell  ij  is  e^j. 

The  pressure  P  is  measured  at  points  'where  several  volume 
elements  meet.  For  our  example,  pressures  ^Ij+l 

Pl+u+l  "*  ttt  respective  points  x^,  y^j  y^;  etc. 

This  is  shown  in  Fig,  3*  Pressures  at  intermediate  points  are 
assumed  to  be  given  through  linear  interpolation  of  nearby  pressures. 
Thus  the  pressure  at  a  point  y  along  face  Oij  is  given  by 


pP  m  p® 

4  *4  4 


l :  yu 

yU+i  ■  y: 


U  "  *ij) 
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A  similar  relation  applies  to  each  of  the  other  faces  lij,  2iJ,  and 

31J. 

External  forces  are  treated  in  much  the  same  manner.  We  observe 
that  the  potential  function  #  appearing  in  Eqa.  (l)  and  (9)  is  always 
in  the  form 


s”4  *  "  (S*  t?  Sz)  ■  (‘x»  -x» -z) 

where  X,  Y,  and  Z  are  the  components  of  force  in  the  x,  yf  and  2. 
directions,  respectively.  For  our  example,  gravity  is  the  potential 
force  and  its  components  are 

grad  «  *  (0,  g,  0)  (19) 

where  t  is  defined  by  Eq.  (12).  In  this  particular  case,  grad  I  is 
a  constant. 

INITIAL  CONDITIONS 

Ohe  mass  in  a  given  element  of  volume  is 


M 


(20) 


If  ve  substitute  Eq.  (l4)  for  p  and  (x  dx  dy  dfl)  for  dV  in  Eq.  (20), 
the  total  watss  M^j  originally  in  cell  ij  is 


Jx11+1 

xli 


/lj+l 

JyU 


Jo  po 


e-y/x 


x  dx  dy  de 


which  after  integration  becomes 
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Mfj  -  ¥  (Xll+1  •  4)  (e'YljA  -  «'WX)  (21) 

This  mass  is  divided  equally  among  the  mass  points  originally  in  the 
cell.  The  points  in  each  cell  are  positioned  so  that  the  density 
is  approximately  uniform  over  the  volume  of  the  cell.  For  instance, 
if  there  are  eight  points  original ly  in  a  given  cell,  their  locations 
might  appear  as  indicated  hy  the  cross  marks  in  Fig.  4.  Since  the 
points  are  free  to  move  from  cell  to  cell  as  time  goes  on,  it  is 
necessary  to  record  the  mass  of  each  point  so  that  the  toted,  mass 
or  density  in  a  cell  can  be  determined  at  any  time. 

The  initial  velocity  in  each  cell  is  zero  in  the  example 

ujj  -  0  (22) 

The  disturbed  area,  within  which  the  specific  internal  energy 
eij  “  a  eQ>  consists  of  a  number  of  cells  which  lie  wholly  or  par¬ 
tially  within  the  radius  rQ.  Outside  this  region,  e^j  *  eQ. 

EQUATION  OF  STATE 

Having  specified  the  other  initial  conditions,  we  calculate  the 
pressures  from  the  equation  of  state,  Eq.  (4).  In  order  to  evaluate 
the  pressure,  which  is  a  function  of  density  and  specific  internal 
energy,  it  is  necessary  to  define  an  effective  density  and  specific 
internal  energy  for  the  region  Immediately  surrounding  the  point  at 
which  pressure  is  computed.  Since  the  density  of  the  fxuid  is  con¬ 
sidered  to  be  constant  over  a  cell,  we  define  the  mass  M  associated 
with  the  pressure  as 


-32- 


-33- 


M  “  (MOiJ  +  ^l-lj  +  M31-1J-1  +  “sij-l1  (23) 

where  the  first  subscript  on  each  M  refers  to  a  subcell  (see  Fig. 
2b).  The  volume  in  which  this  mass  is  contained  is  given  by 

V  ■  2  (x3i  '  x3i-l)  (y3J  ”  y3J-l)  ^ 

Substituting  in  Eq.  (17),  we  find  that  the  density  associated  vith 
the  pressure  is 


(M< 


'oij  *  *  M3i-1.1-l  *  **2: 

l(X31  '  x3i-l)  (y3J  "  y3J-l 


W) 

) 


(25) 


In  a  similar  manner,  the  effective  specific  energy  for  the 
pressure  computation  becomes 


M01.1e1.1  +  “ll-lA-lJ  +  M21J-lelJ»l  *  M3i-lJ-lel-U-l  (26) 


Upon  substitution  of  Eqs.  (25 )  and  (26)  in  the  equation  of  state 
for  the  fluid,  Eq.  (ll),  the  pressure  is  determined. 

Certain  boundary  conditions  affect  the  computation  of  pressure. 
For  instance,  at  y^j  a  0,  either  Eqs.  (23)  through  (26)  must  be 
modified  to  take  into  account  the  presence  of  the  boundary,  or  mass 
and  specific  energy  for  fictitious  cells  below  y  =  0  must  be  supplied. 
Similarly,  some  adjustment  is  needed  along  the  axis  of  symmetry  at 
xli  "  °*  We  choos^  to  modify  the  equations  rather  than  supply 
variables  for  the  fictitious  cells.  Hence,  at  «  0,  y ^  -  0, 

Eqs*  (23)  through  (26)  become 
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M  -  M^j 

V  .  I  (4)  (»M) 

p  "  1  s 

I  x31y3J 

e  - 

vhere  far  this  example  -  0,  -  0  vhen  i  -  1  and  j  -  1 


*U  *  °»  yu 

>  0,  Eqs.  (23)  through  (26)  become 

M  "  (M01J  +  “ll-lj) 

(23b) 

V  - 1(4  -  4-i)  fo) 

(24b) 

(M01.1  +  “u-u) 

1  (4  ■  4-0  (yij) 

(25b) 

-  -  M0ifl.1  +  "li-lA-M 

M01J  +  “li-u 

(26b) 

Finally,  at 

=,  C,  y1,  /  o,  Eqs.  (23)  through  (2b)  become 

«  ■  ("on  *  "21J-O 

(23c) 

»  ’  1  («3»)  (y3J  '  y3J-l) 

(24c) 

(23a) 

(24a) 

(25a) 

(26a) 
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p  -  - Wgl^'1)-  (29o) 

2  X3i(y3J  “  y3J-l) 

e  .  M01Jei.J  +  (a6c) 

M01J  eij-l 


Other  types  of  boundary  conditions  require  similar  considerations. 
Free  boundaries  are  such  that  the  pressure  must  be  zero  on  the  bound¬ 
ary.  This  can  be  handled  numerically  by  imposing  the  following 
restriction:  if  any  is  zero  for  the  full  cell  corresponding  to 
in  Eqs.  (23),  (23a),  (23b),  or  (23c),  whichever  is  applicable, 
then  the  resulting  pressure  is  zero;  otherwise  the  computation  of 
pressure  is  carried  out  as  previously  described. 

Ohls  method  of  evaluating  pressure  from  the  equation  of  state 
is  employed  not  only  in  the  computation  of  initial  conditions  but 
also  in  subsequent  computation  in  the  integration  process  itself. 


PRELIMINARY  CALCULATION  OF  VELOCITY 

The  first  step  in  the  numerical  integration  of  the  equations  is 
to  estimate  the  new  velocity  in  each  cell.  The  x- component  of  Eq. 
(l)  is  written 


Du 


+ 


&P 

5x 


c 


which,  for  our  example  is 


Du 


+ 


3P 

3x 


0 
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since  grad  *  in  the  x  direction  is  zero.  Neglecting  the  convective 
tern  and  dividing  by  p  we  have 


bu  m  1  &P 

EE  "  p  3x 


(27) 


>p 

This  equation  requires  p  and  to  be  known  in  order  to  evaluate  the 

change  in  the  u  component  of  velocity  with  respect  to  tine.  For 

dP 

the  fluid  model  the  density  p  is  known  in  each  cell,  but  ^  is  not. 

xp 

We  estimate  ^  in  the  following  manner. 

Assuming  that  the  pressure  varies  linearly  from  one  corner  of 
a  cell  to  another,  an  average  pressure  P  for  a  face  can  be  computed. 
Far  f&ce  0  of  cell  ij  the  average  pressure  P^j  is 


(  TS 


(28) 


and  for  face  1  the  average  pressure  P^j  is 


^iij  “  i  (^u  +  ^ij+i ) 


dp 


If  we  approximate  with 


(29) 


^  *11+1  "  Xli 


and  substitute  this  expression  into  Eq.  (27),  we  have  far  the  change 
in  u  in  cell  i S  at  time  t 

n 


-i  :  *0^. 
p^v^n+i  -  X11/ 


(30) 
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vhere  the  dot  denotes  the  tine  derivative  n»e  pressures  and 
are  given  hy  Eqs.  (28)  and  (29)  respectively,  and  p°j  is  given 
by  Eq.  (17). 

For  the  change  in  velocity  in  the  y  direction,  the  differential 
equation  is 


Dv  ,  3p  .  ,  .  ,3* 

*Di  +  ^+  1  + 


which  for  our  example  may  be  written 


^  +  S  + « - 0 


Again,  vc*  drop  the  convective  term  and  divide  by  p,  to  obtain 


dv  -1  dP  „ 

St  ■  pSy  '  s 


(31) 


Calculating  average  pressures  in  the  manner  stated  above,  we  have 
for  cell  ij 


1J 


'U*!  -  } 


6 


where 


f2iJ  ■  5  (^J  +  *1+13  ) 


p“  ,±/y  +  p“ 

3iJ  2  V  U+1  i+lj+l 


) 


(32) 


(33) 

(34) 


If  the  change  in  velocity  is  assumed  to  be  constant  over  the  time 
step  and  equal  to  the  value  at  the  beginning  of  the  time  step, 
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tlw  «>tia»t«d  velocities  at  are 


U1J  “  U1J  +  ^n+l  *1J 


^lj1  “  V1J  +Afco+l  *lj 


(35) 


with  m  -  t^.  Here  the  tilde  denotes  an  estimate  of 

velocity,  since  we  have  astuaed  that  masses  do  nut  move.  Also,  the 
estimated  average  velocities  for  the  time  step  are 


-n+1 

1 

/  n 

.  *wl' 

Uij 

"2 

Vuu +  uij  , 

-n+1 

■i 

f  n 

.  ^nfl 

ViJ 

(36) 


H1EUMIHARY  CAUCULATIOH  OF  DITHUUL  jgggOY 

The  next  phase  of  the  numerical  process  la  concerned  with  the 
estimation  of  specific  Internal  energy  in  each  cell  at  the  end  of 
the  time  interval.  The  conservation- of-energy  relation  given  by 
Ed*  (9)  serves  as  the  basis  for  this  computation.  Since  we  assume 
that  density,  velocity,  and  specific  internal  energy  are  constant 
over  the  cell,  and  grad  t  Is  also  constant,  Eq,  (9)  may  be  written 


pV 


u  •  grad  •  u  +  Pu  •  dS 


Dropping  the  convective  term  end  dividing  by  pV,  we  have  for  the 
change  in  specific  energy 

(37) 
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Integratlmg  Eq,  (37)  with  respeot  to  time  yield*  the  change  In  spe¬ 
cific  energy  for  the  time  Interval,  which,  when  added  to  the  specific 
energy  at  the  beginning  of  the  interval,  results  In  the  desired  esti¬ 
mate  of  specific  energy*  Because  of  this  Integration  with  respect 
to  time,  average  velocity  for  the  time  step  Is  used  wherever  velocity 
Is  required  In  evaluating  Eq.  (37) • 

For  the  physical  example  vhere  grad  t  Is  given  by  Eq,  (19),  the 
first  tens  In  Eq.  (37)  Is 

-u  •  grad  *  -  -gvJJ1  (38) 

and  the  second  tens  Is 

-jt(5*-5)  (39) 

The  mass  pV  In  the  third  tern  In  Eq.  (^6)  Is  equal  to  V^,  while 
the  Integral  represents  the  rate  at  vfaich  the  energy  Is  btli^  dimin- 
Ished  by  work  done  on  the  surface  of  the  volume  element.  In  evalu¬ 
ating  the  surface  work  integral  for  a  cell  we  note  that  the  voluee 
element  under  consideration  has  six  different  surfaces.  Only  four 
of  the^e,  however,  contribute  to  the  integral  since  the  net  work 
done  on  the  surfaces  coincident  with  the  meridian  planes  Is  zero  due 
to  cylindrical  symmetry.  Integrals  for  the  four  other  surfaces 
(0,  1,  2,  3)  which  do  contribute  are  evaluated  separately  then  susssed 
to  determine  the  total  integral. 

To  evaluate  each  of  these  Integrals,  we  first  determine  the 
quantity  Pa  *  d S  and  then  Integrate.  For  face  01 J  the  pressure  P 


-39- 


Integrating  Eq.  (37)  vith  respect  to  tine  yields  the  change  In  spe¬ 
cific  energy  for  the  tine  Interval,  vhloh,  when  added  to  the  specific 
energy  at  the  beginning  of  the  interval,  result*  in  the  desired  estl- 
nate  of  specific  energy.  Because  of  this  integration  vtth  respect 
to  tine,  average  velocity  for  iiie  tine  step  is  used  vherevsr  velocity 
is  required  in  eval  utiag  Eq.  (37)* 

Far  the  pi  ic&l  example  vhere  grad  #  Is  given  by  Eq.  (19),  the 
first  term  in  Eq.  (37)  is 

-u  •  ©tad  *  -  -gv^1  (36) 

and  the  sec©*-?'  tens  la 


.5^1  ..  v^1  V1 

wij  u  vij  ^ 


<n) 


The  mu  pV  In  *'  talrd  ten  In  R*.  (36)  Is  equeu.  to  p“^  411.  s 

the  integral  represents  the  rate  at  uhlch  the  energy  is  being  dimin¬ 
ished  by  work  done  on  the  sw  tace  of  the  volxarc  element.  In  evalu¬ 
ating  the  surf1  *  work  integral  for  a  cell  ve  note  that  the  volimtt 
element  undex  consideration  has  six  different  surfaces.  Only  four 
of  these  hr* eve/,  eontrluute  to  the  integral  since  the  net  work 
done  the  surfaces  c  incident  vith  the  meridian  planes  is  zero  due 
wo  cylindrical  syi r  -ry.  Integrals  for  the  four  other  surfaces 
(0,  1;  c,  3)  vhicV.  do  contribute  are  evaluated  separately  then  sunned 
to  determine  the  total  integral. 

To  evaluate  each  of  these  integral s,  ve  first  determine  the 
quantity  Pu  •  ^3  and  then  integrate.  Far  face  PiJ  the  pressure  P 


is  a  function  of  y  alone  as  given  in  Bq,  (16)  •  The  function  u  *  df} 
represents  the  velocity  nurwal  to  the  surface  (and  taken  to  he  posi¬ 
tive  in  the  outward  sense)  tines  the  ulesnnt  of  arva  over  which  It 
acts.  Since  ve  consider  the  velocity  in  each  cell  to  he  ooastast 
over  the  cell,  a  discontinuous  velocity  at  a  surface  emewr  to  two 
cells  can  ocjut  If  the  velocities  in  those  cells  are  different.  Hence 
the  following  assumption  Is  eade:  Die  velocity  on  such  a  coion 
surface  is  equal  to  the  average  of  the  velocities  in  the  adjacent 
cells.  The  function  u  •  d3  nay  now  be  evaluated;  for  face  0  It  is 

u  •  as  -  -  |  (%  ♦  tiff)  **)  (ko) 


Die  negative  sign  appears  because  of  the  outward  sense  of  the  noraal. 
Upon  substitutions  of  Eqs.  (18)  and  (to)  ve  find  for  this  one  face 


which  after  integration  bee  ones 

Pu  •  dS  -  ^jp  (yxj+x  -  Yx +  +  “ij4")  ^ 

Proceeding  in  a  similar  wanner  for  the  opposite  face,  lij,  we  find 
that  the  integral  of  Eq.  (37)  ia 

- J  PU  •  dS  -  -  ^i(yli+1  -  *  I1+1J+1)(UJ1  + 


Next,  the  pressure  along  face  21J  Is  given  by 


\  (*»«  •  *«)  (W) 

and  the  function  u  *  dS  in  the  outward  sense  is 

-u  •d3--i(v^1  +  ^J1)xdx  (>*) 

Substituting  for  the  integral,  we  have 


which  after  integration  becomes 

^  ^  “  +  15  (*u.l  ‘  xll)  +  vl'j1)  j*lj  (^ll  +  *li+l) 

*  *?.«  (xu  *  aw)]  <«> 

In  a  similar  fashion  for  the  integral  for  the  opposite  face,  3iJ,  is 
pS  •  dS  -  -  35  (xn+i  -  xu)  (^j1  +  vlJ+l)  [^j+l  (^li  +  xn+i 

+  +  2xU+l.)j  ^ 

The  total  integral  in  Eq.  (37)  is  the  sum  of  the  integrals  for  the 
four  faces  given  by  Eqs,  (4l),  (42),  (4j),  and  (46U 

It  should  be  roted  that  when  one  cell  receives  an  increment  of 

H*  -* 

energy,  Pu  •  dS  At,  exactly  that  amount  is  subtracted  from  the 


appropriate  adjoining  cell,  so  that  total  energy  is  conserved.  Since 
the  re  partitioning  process  discussed  in  the  following  sections  also 
conserves  energy,  the  external  potential,  #,  is  the  only  agency  %hich 
can  cause  a  change  in  the  total  energy.  In  cases  where  ♦  *  0,  this 
provides  a  waivable  check  on  the  computation. 

The  change  In  specific  energy  per  unit  tine,  Eq.  (37),  may  now 
be  evaluated  using  Bqs.  (3B),  (39),  (4l),  (42),  (45),  (46)  *ith 

pV  *  Integrating  with  respect  to  tine  we  find 


Htt^l 

etj 


eU  ' 


n+1  A. 
^4  -*■ 


nH 


Os 


-n+l.n  ^  -n^l.n  \ 
“14  +  14  14/ 


*1+14+1 1*!!  *  2xli+l})  "  ^vi4-l  +  *14*)  ” 

xu+l}+  P1+14{X11  +  2xU+l}) 

l 

In  Eq.  (by)  no  account  has  been  taken  of  the  effect  of  boundary  condi¬ 
tions  other  than  the  presence  of  free  boundaries  where  the  pressure 
is  zero  (as  discussed  previously).  In  our  particular  example  no  work 
can  be  done  on  the  fixed  boundary  or  on  the  boundary  of  s>f  try 
because  the  velocity  normal  to  these  boundaries  Is  zero.  Hence  the 
ter*  in  Eq.  (47)  resulting  from  Eq.  (45)  is  zero  far  y^  » 
the  term  resulting  fro*  Eq.  (4l)  is  zero  when  -  0. 


0,  and 


Having  computed  the  estimates  for  the  velocity  and  specific 


internal  energy  at  time  t  ve  can  determine  the  momentum  in  each 
direction  and  the  sura  of  the  kinetic  and  interned  energy  for  each 
cell*  The  momentum  in  the  x  and  y  directions,  respectively,  are 


Momeatumx  « 

Momentum^  -  M^j  v^1 


m 


where 


<3  ■  »?3 


The  estimated  kinetic  plus  internal  energy  in  a  cell  is 


Expressions  (48)  and  (49)  are  correct  at  the  end  of  the  time  interval, 
provided  no  mass  changes  location  from  one  cell  to  another  (i.e., 

■c  "  ^ij)  *  correc*  these  estimates  of  momentum  and  energy  for 
cells  which  do  not  meet  this  condit-<™\  each  mass  which  changes  cells 
is  considered  to  carry  with  it  to  th*  new  cell  its  share  of  the 
momentum  and  kinetic  plus  internal  energy.  This  together  with  the 
integration  of  Eq.  (37)  insures  that  mass,  momentum,  and  toted  energy 
(kinetic  +  internal  +  potential)  are  all  conserved. 


MASS  MOVEMENT 

In  the  next  phase  of  the  computation,  the  masses  are  moved 
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ac cording  to  an  average  velocity  for  the  time  cycle.  First  the  sub¬ 
scripts  p  and  q  are  determined  such  that 


vhere  x and  y11  define  the  position  of  the  mass  point  at  time 
m  m 

tn  and  and  y2^  define  the  center  of  area  of  the  cell  pq.  Hie 
velocity  of  the  mass  thus  located  is  considered  to  be  a  sum  of  the 
velocities  in  cells  p-1  q-1,  p  q-1,  p-1  q,  and  p  q  weighted,  respec¬ 
tively,  by  aQ,  ax,  ag,  and  a3 


vhere  and  vjj*1  are  the  migration  velocities  assigned  to  the 
mass  point  in  the  x  and  y  directions,  respectively.  The  weights  are 
determined  in  the  following  fashion:  First  the  proportional  areas 
Cq,  qv  c2,  and  c^  are  computed 
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Flnally,  the  weights  aQ,  a^,  a,  ,  and  a^  are  defined 


a0  *  t  +  bib  +  b 
0  12  3 


1  b0  +  bl  +  b2  +  b3 


(5»0 


“2  -  t0  VT-+  b2V-b3 


a  ,  3 _ 

3  +  ^1  +  ^2  +  ^3 

Substituting  Eqs.  (54)  in  Eqs.  (51)  yields  the  average  velocity  far 
the  tine  cycle  far  mass  m. 

Again,  boundary  conditions  modify  the  computation  of  velocity. 
For  the  fixed  boundary  at  y  -  0,  cells  p-X  q-1  and  p  q-1  represent 
fictitious  cells.  Die  determination  of  the  velocity  near  this 
boundary  is  achieved  by  modifying  the  computation  of  the  c»s  of 
Eq».  (52) 


co 


0 


C 


1 


0 


(52a) 
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Substitution  of  Eqs.  (52a)  in  Eqs.  (52)  and  continuing  in  the 
manner  previously  described  permit  the  velocity  of  the  mass  near  the 
boundary  to  be  determined. 

Similarly,  for  the  velocity  near  the  boundary  of  symmetry,  the 
c 1  s  are  computed  in  the  following  manner. 


c0-° 


K  V  *20  -  ym  \ 

1  \x2pA 


Cg  “  0 


(52b) 


'3  U 


'2p  I 


ym  -  y2q-l 

^2q  ■  y2q-l  / 


For  the  velocity  near  x  «  0  and  y  *  0  the  c*s  are 


c0  =  ° 


c 


1 


0 


°2  =  0 


(52c) 


c 


3 


The  new  location  of  the  mass,  given  by  and  is 

determined  by  integrating  Eqs.  (51)  with  respect  to  time 


n+X 


n  -rH-1 

x  *  A5:  1  u 

■  n*l  o 


+  at 


n*l  n 


r^l 


(55) 


Furthermore,  we  Impose  boundary  conditions  at  x  =  0  aod  a;  y  *  0  such 

x 

that  no  Bass  can  cross  the  boundary.  Thus  If  x  <  0,  x  is  set 

n  si 

equal  to  fxf*1 1  similarly.  If  y**1  <  0,  yj**1  is  set  equal  to  far***1 1. 

D  n  D  B 

The  n*1  mss  point  at  tine  t  In  In  cell  l°jJ,  itee  1®  and  J°  are 
defined  as 


‘I-1  **“  *11  *  x»  <  *11*1 

(56) 

If  at  tine  t^,  and  determined  by  substituting  x£**  and 

y15*1  for  x»  and  yn  in  Eq.  (56),  are  equal  to  i®  and  j“  respectively, 

9  B  B  B  3r 

no  adjustaent  in  Boaeotu  and  energy  is  required  since  the  nass  does 

not  change  cells.  If,  however,  the  mass  does  change  cells,  then 

estimates  of  aonentts  and  energy  In  both  the  old  cell  and  the  nev 

th 

cell  are  aodlfied.  The  estimated  aaomtua  of  the  n  mss  point 

for  the  tiae  interval  t  ^  t  s  t.  . ,  is,  in  the  x  direction  and.  7 
el  mi 

direction,  respectively 


■tfS1- 


M  V^1 

Vu 


vhere  Is  the  nass  of  the  a^  nass  point,  the  estl—ted  total 
energy  assigned  to  this  saae  aass  point  Is 


(57) 
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M11 

r*1*  m 

^ij  2 


+  7*1 


s1*)} 


(56) 


where  e^J\  The  total  mass  in  the  cell  at  time  tfi 

is  is  the  total  estimated  internal  energy  fo;:  this 

same  cell  before  mass  movement.  Notice  that  i  =  i*J  and  )  =*  in 

m  m 

both  Eqs.  (57)  and  (58).  In  going  from  cell  i*j£  to  the 

mass  takes  with  it  both  momentum  and  energy.  Hence  expressions 

(57)  and  (58)  are  subtracted  from  the  estimated  momentum  and  kinetic 

plus  internal  energy  far  cell  i^J^  and  added  to  those  of  cell 

for  each  mass  point  which  crosses  from  cell  i*J**  to  cell 

r  m  m  m  m 


REPARTITIONING 

After  ell  of  the  masses  have  been  moved  and  the  total  estimated 
momentum  ar*i  energy  for  the  cells  have  been  adjusted  correspondingly, 
a  final  set  of  values  for  density,  velocity,  and  specific  internal 
energy  at  time  ntl  must  be  assigned  to  each  cell.  These  values  are 
chosen  so  that  the  total  mass,  momentum,  and  energy  assigned  to  each 
cell  after  mass  movement  are  conserved. 

First,  the  mass  of  the  ijth  cell  after  mass  movement  is 

found  by  simply  adding  up  the  masses  of  all  of  the  individual  mass 
points  in  the  ij^  cell  at  this  phase  of  the  calculation,  lhe  nev 
density  of  the  cell  is  then  given  by 


Mn+1 

„n+i  =  £l 


ij 


(59a) 
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The  new  velocities 


/  n+1  n*-l\ 

V*1J  '  viJ  ) 


are  given  respectively  by 


[Total  estimated  momentum,  x  direction]^ 


1J 


w^+1 

"Ij 


(59b) 


[Total  estimated  moment,®,  y  direction], , 

•ft1 - pi - u  <»» 

"ij 


nfl  [Total  «timated  energy]^  -  |  mJJ1  [(u^1)2  +  (v^1)2] 

e  - 

“ij 


(59d) 


mESSURE  CALCULATION 

The  final  step  in  the  numerical  Integration  process  Is  to  compute 
the  pressures  as  described  in  Section  IV  using  the  new  values  of 

specific  internal  energy  e^1  and  density  p^j1. 

STABILITY  CHECK 

Together  with  the  computation  of  the  internal  pressures  it  is 
desirable  to  perform  a  check  on  the  stability  of  the  numerical  inte¬ 
gration  process.  At  the  present  time,  th'ire  appears  to  be  no  sta¬ 
bility  criterion  specifically  designed  for  the  PIC  process,  although 
the  von  Neumann  method  ^  of  stability  analysis  used  in  associa¬ 

tion  with  the  existing  PIC  analogue  difference  equations^'®)  could 
probably  be  used  to  develop  one.  A  rough  approximation  to  a  sta¬ 
bility  criterion  is  given  by  the  conditions  that  the  mesh  speed 
in  any  direction  should  exceed  the  speed  of  a  sound  wave  (relative 
to  a  fixed  reference  frame)  in  that  direction.  For  our  problem  this 
condition  becomes, 


Min  i  _  —Ay  \ 

l|u|+  V  lVi+  VS  J 


where  Ax  and  Ay  ^re  increments  of  Euler i an  coordinates  and  v  is  the 

s 

local  velocity  of  sound  of  an  observer  moving  with  the  fluid. 

This  coalition  is  always  met  by 


At  <  Min 


v  +  “  2 

s  VU  +  V 


For  practical  c amputations  we  have  found  that 


/u  *►  v2  s  v 


We  therefore  choose  as  our  stability  criterion 


| ax,  Ay 


where  r  is  a  constant  ( approximately  equal  to  2)  that  is  generally 
established  by  preliminary  computational  experiments  for  whatever 
problem  ve  have  in  mind. 

The  velocity  of  sound  is  given  by 


/  dp\ 

=ws 


that  is,  the  partial  derivative  of  P  with  respect  to  p  at  constant 
entropy.  Sinco  at  constant  entropy 


de  =  -§  dp 
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then 


2 


v 


3 


(62) 


For  stability  of  the  numerical  process  to  be  assured,  the  C  cur  ant 
condition,  Eq,  (60),  is  observed  thr cutout  the  region.  In  the  case 
of  the  equation  of  state  given  by  Eq*  (10)  we  have 


2 


V 


3 


(63) 


Therefore  for  each  point  in  the  fluid  at  which  pressure  is  computed 
the  condition  to  be  satisfied  is 


s  ■  r 


(64) 


In  addition  to  a  check  on  stability  we  may  use  Eq.  (64)  to  adjust  the 
time  interval  during  the  course  of  integration  so  that  an  optimum 
time  interval  is  always  used.  This  is  done  in  the  following  manner 


"*•2  =  Atnfl  f0r  S1  s  s  s  s2  <  1 


-  riAVi  for  s  <  S1  (ri  >  x) 

AVa  =  far  s2  <  3  (r2  <  X) 


That  is,  if  the  calculated  value  of  the  stability  function  s  at  time 

t* ,  lies  between  or  on  either  of  two  preassigned  numbers,  s  <  sn  <  1, 
n+l  2 

then  the  subsequent  time  interval  is  not  changed.  If  the  stability 
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function  s  is  less  than  s^,  then  the  time  interval  is  increased  by  a 
multiplier  Ty  If  p  is  greater  than  s2,  hoover,  the  tine  interval 
is  decreased  by  the  multiplier  r^.  The  numbers  r^  and  rg  are  also 
positive  pre as signed  constants, 

GRID  CHARGE 

Sooner  or  later  the  disturbance  induced  by  the  presence  of  the 
blast  vill  reach  the  grid  boundary  at  x  «  xmax>  7  =*  or  both. 

At  this  point  in  time,  computation  with  the  current  model  must  cease, 
because  the  material  will  soon  flow  to  nonexistent  cells  and  further 
computations  will  become  meaningless.  Consequently,  at  the  end  of 
each  integration  step,  tests  are  made  to  determine  whether  or  not  the 
disturbance  has  reached  the  boundary.  These  tests  involve  a  search 
for  nonzero  velocity  and  far  specific  internal  energy,  density,  or 
pressures  which  are  different  from  their  original  values  in  any  cell 
adjacent  to  the  boundary. 

In  instances  where  disturbances  spread  through  a  larger  and 
larger  volume  of  material,  it  is  often  desirable  to  maintain  at  all 
times  a  high  degree  of  resolution  with  respect  to  the  disturbed 
volume.  The  airburst  problem  is  one  such  case.  At  the  time  when 
a  disturbance  reaches  the  assumed  boundary  in  the  example,  it  is 
possible  to  create  a  new  model  with  the  characteristics  of  the  old 
model  at  that  specific  time,  so  that  integration  can  be  resumed. 

This  is  done  in  the  following  manner. 

The  new  region  is  selected  in  such  a  way  that  all  the  fluid  in 
the  old  model  is  contained  within  the  volume  of  the  new  model.  Such 
selection  for  our  physical  example  is  illustrated  in  Fig.  5*  Here 
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the  maximum  dimensions  of  the  nev  region  are  at  "new”  x  and  "new" 

max 

ym  :  the  old  region  is  completely  enclosed  in  the  nev  region, 

UlfcLX 

The  new  region  is  then  subdivided  into  volume  elements  in  the 
same  manner  as  described  before.  Mass,  momentum,  and  energy  are  now 
assigned  to  the  cells  of  the  whole  nev  region  in  two  distinct  steps: 
First  to  those  cells  or  parts  of  cells  which  coincide  with  the  nev 
region,  and  then  to  those  cells  or  parts  of  cells  which  coincide  with 
the  old  region.  Initial  values  of  mass,  kinetic  plus  internal  energy, 
and  momentum  are  provided  in  the  new  cells  for  the  additional  fluid 
which  has  been  encompassed;  i.e.,  in  the  new  cells  where  >  "old" 

xmax'  030  :flj+l  >  ymax*  Here  is  P°ssitlje  'that  a  given  cell 

may  encompass  some  old  volume  as  veil  as  some  nev  volume  of  fluid. 

In  this  event,  the  mass,  energy,  and  momentum  for  the  new  fluid  only 
are  provided  in  these  cells. 

Next  the  mass  of  the  fluid  in  each  cell  of  the  old  region  is 
added  to  the  proper  nev  cells,  as  shown  in  Fig.  6,  under  the  assump¬ 
tion  of  uniform  density  for  each  old  cell.  The  portion  of  mass  of 
the  old  cell  which  is  contained  in  the  new  cell  1  is  added  to  the 
mass  in  new  cell  1.  Similarly  the  kinetic  plus  internal  energy  and 
momentum  associated  with  mass  are  added  to  the  energy  anfl  momentum 
in  nev  cell  1.  The  other  nev  cells,  2,  3,  and  4,  receive  mass  M^,  M^, 
and  respectively,  from  the  old  cells  and  Oso  receive  the  corre¬ 
sponding  portions  of  energy  and  momentum. 

After  all  the  mass,  energy,  and  momentum  in  the  old  cells  have 
been  redistributed  to  the  nev  cells,  the  velocity,  specific  energy, 
and  density  are  computed  in  the  new  cells  in  tne  manner  described  in 
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Section  IV.  The  mass  in  each  cell  is  then  divided  equally  among  a 
number  of  mass  points,  positioned  so  that  the  density  Is  approximately 
uniform. 

Finally  the  pressure  at  the  corners  of  the  cell  is  computed  from 
the  equation  of  state  of  the  fluid  using  the  method  detailed  in  Sec¬ 
tion  IV. 

The  variables  associated  vith  the  nev  cells  are  the  Initial 
conditions  for  continuing  the  integration.  As  time  goes  on,  subse¬ 
quent  grid  changes  may  be  necessary  to  obtair  the  desired  solution 
to  the  physical  problem. 

The  grid  change  pro /Ides  a  means  by  which  we  can  overlay  on  a 
given  network  one  of  less  detail,  while  conserving  mass,  momentum, 
and  energy.  Some  Joss  of  detail  occurs,  to  be  sure,  but  the  general 
characteristics  of  the  old  region  arc  carried  over  into  the  nev  region. 
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